Hierarchical power control of a large-scale wind farm by using a data-driven optimization method

With the participation in automatic generation control (AGC), a large-scale wind farm should distribute the real-time AGC signal to numerous wind turbines (WTs). This easily leads to an expensive computation for a high-quality dispatch scheme, especially considering the wake effect among WTs. To address this problem, a hierarchical power control (HPC) is constructed based on the geographical layout and electrical connection of all the WTs. Firstly, the real-time AGC signal of the whole wind farm is distributed to multiple decoupled groups in proportion of their regulation capacities. Secondly, the AGC signal of each group is distributed to multiple WTs via the data-driven surrogate-assisted optimization, which can dramatically reduce the computation time with a small number of time-consuming objective evaluations. Besides, a high-quality dispatch scheme can be acquired by the efficient local search based on the dynamic surrogate. The effectiveness of the proposed technique is thoroughly verified with different AGC signals under different wind speeds and directions.


Introduction
In respond to the energy shortage and environment degradation, the renewable energy like solar energy [1] has gained a rapid development in recent years around the world.As the backbone of renewable energy, the global new installed capacity of wind energy has grown remarkably, where the year-over-year growth is more than 50% from 2019 to 2020; and the compound annual growth rate of the offshore wind energy will be more than 30% from 2021 to 2025 [2].In general, the power output of renewable energy is highly influenced by the weather conditions, such as wind speed, solar irradiance and temperature [3].With the rapid increasing renewable energy, the power system is facing tremendous pressure on the real-time power balance due to the decreased moment of inertia and the increased random power fluctuation [4].Hence, it requires more and faster reserve resources to participate in automatic generation control (AGC).
Since a large-scale wind farm (WF) can provide a sufficient regulation capacity and a superior regulation performance, extensive studies on AGC have been investigated by considering the deep participation with WF.To improve the dynamic regulation performance of the control area, a lifelong learning [5] was proposed for coordinated control between WFs and other AGC units.In [6], an adaptive distributed auction-based algorithm was presented to achieve a distributed optimization of AGC dispatch.Based on the similar control framework, a deep reinforcement learning [7] was designed to coordinate WFs and other energy resources in an integrated energy system.All of these works mainly focused on the control process of the whole system, in which a WF was regarded as a single wind turbine (WT).In fact, a WF consists of numerous WTs, where the AGC signal of the whole WF should be distributed to each WT.
From the viewpoint of WF, a set-point tracking control [8] was designed for AGC, in which the power command of each WT is proportional to its available power capacity.This simple proportional dispatch method can not satisfy the optimal control requirement of WF for AGC.To handle this issue, a distributed low-complexity controller based set-points optimization [9] was presented for fatigue reduction of WF.In [10], a distributed model predictive control was proposed for minimization of the WT load and improving the tracking performance of WF.Besides, an alternating direction method of multipliers (ADMM) based hierarchical optimal power control [11] was presented to guarantee the rotor speed stability, minimize the energy loss, and damp the oscillations.In [12], the distributed reactive power control can be achieved for the large-scale WF cluster by ADMM.Like the proportional dispatch, a leader-follower consensus protocol based on the fair load-sharing [13] was employed to realize a distributed control of power output regulation for WF.Similarly, a multi-agent framework based cooperative output regulation [14] was designed for a large-scale WF following leader-follower consensus control.Moreover, an efficient distributed economic model predictive control strategy [15] was proposed to take the power tracking and economic operation into account for a large-scale WF.To effectively reduce the optimization complexity and difficulty, the clustering algorithm [16] is usually used for dividing the large amount of WTs.In [17], an incremental clustering based equivalent method was proposed to improve the equivalent accuracy for the large-scale WFs.In [18], the large-scale WFs were classified into multiple groups by using a dynamic strategy, thus the real-time dispatch can be easily solved based on the grouping results.Based on the WTs clustering, a hierarchical power dispatch [19] was designed for a large-scale WF based on the proportional control and proportional distribution strategy.In [20], the hierarchical power control was extended to multiple large-scale WFs with active and reactive power control concurrently.However, these works did not consider the interactive influence between different WTs by the wake effect [21], i.e., the inlet wind speed of the downstream WTs can be influenced by the operating state of the upstream WTs.
By taking the wake effect into account, various methods [22] were proposed for power control of WF, including particle swarm optimization (PSO) [23], Bayesian Ascent based datadriven optimization [24], model-predictive optimization [25], cooperation-driven distributed control [26,27], double-layer machine learning [28], and so on.Most of these studies essentially belong to a single-degree-of-freedom control based dispatch, where the optimization space is restricted since the rotor speed control (RSC) and pitch angle control (PAC) are separated.To simultaneously implement these two control degrees, a small-signal method based two-degree-of-freedom control was constructed in [29], which can reduce the fatigue load during the power tracking.Following the same manner, a clustering-based coordinated control method [30] was proposed for power system frequency support by combining the kinetic energy, but it did not address the power tracking via distributing the AGC signal.
On the whole, the above studies did not consider the dynamic regulation performance in the control objective, i.e., a WF should trace the AGC signal as quickly and accurately as possible.In general, the response speed of RSC is faster than that of PAC [31], but the regulation range of WT power output by RSC is much narrower than that by PAC.Consequently, the control commands of these two degrees for all the WTs should be optimized to acquire a satisfactory dynamic regulation performance under different AGC signals.Meanwhile, it easily leads to a high optimization difficulty and an expensive computation due to the large number of controllable variables.To resolve these problems, this paper constructs a new hierarchical power control (HPC) for AGC signal tracking of WF.Like the conventional power dispatch optimization in a WF [23], the control command optimization of each group in HPC is a non-linear and non-convex problem.Compared with the gradient-based algorithms (e.g., the interior point method), the gradient-free algorithms are more competitive and flexible for this optimization as they are highly independent on the specific model of HPC with wake effect.As the most common gradient-free algorithms, genetic algorithm (GA) [32] and PSO [33] have been successfully applied to the optimal power dispatch in WF due to their high global searching ability and easy implementation.Naturally, it can be solved by other meta-heuristic based gradient-free algorithms.But the meta-heuristic algorithms easily result in a long computation time for a high-quality optimum because they should require adequate fitness evaluations.Hence, they are hard to achieve the real-time optimization of HPC with control time cycle from several seconds to minutes, especially for the large-scale WF with numerous WTs by considering the wake effect.
In recent years, the surrogate-assisted optimization [34] has become a popular way to avoid a high computation cost resulted from the large amount of evaluations on the objective or constraint functions.It usually approximates the objective or constraints functions by machine learning models based on the small amount of evaluation data [35], thus a fast optimization can be implemented by different algorithms based on the surrogate model.Inspired by this advantage, a data-driven surrogate-assisted optimization (DDSO) is employed for HPC.Compared with the existing studies, this work has two novelties, as follows: • For the power control of a large-scale WF, most of the existing studies were carried out based on a single-degree-of-freedom control based dispatch, where the optimization space is restricted since RSC and PAC are separated.In [29,30], the constructed two-degree-of-freedom control can simultaneously implement these two control degrees, but they did not consider the regulation performance difference between RSC and PAC.As a result, it easily result in an unsatisfactory regulation performance for AGC.In contrast, the regulation performance difference between RSC and PAC is fully consider in this work, which can dramatically improve the dynamic regulation performance of the whole WF for AGC via twodegree-of-freedom control of rotor speed and pitch angle.
• Compared with the direct power control framework [28], the constructed HPC can effectively the optimization complexity and dimension due to the negligible wake interaction between different WT groups.Moreover, the presented DDSO can result in a much low computation cost against to the meta-heuristic algorithms [23] since it can reduces the number of fitness function evaluations.As a result, a high-quality AGC dispatch scheme can be rapidly acquired for a large-scale WF.
The rest of this work is organized as follows: Section 2 introduces the control framework and optimization model of HPC.The optimization principle and design by DDSO are provided in Section 3. Section 4 gives the simulation results and analysis, and Section 5 concludes this work.

Control framework
To achieve a dynamic power balance between supply and demand, the independent system operator (ISO) should implement the AGC for the control area based on the real-time deviations of frequency Δf and tie-line exchange power ΔP T with other control areas [36], as shown in Fig 1 .In each control interval, the controller of AGC will generate the total dispatch command, then ISO will distribute it to different reserve resources.In respond to the AGC signal from ISO, the control center of a WF will carry out two operations.Firstly, the AGC signal of the whole WF will be distributed to different groups in proportion of their regulation capacities.Then the AGC signal of each group will be balanced by the corresponding WTs via searching the optimal rotor speeds and pitch angles.

WT model
Based on the aerodynamic rule, the energy conversion of the ith WT under different wind speeds can be described as where P i and F T,i are the mechanical power and the thrust force of the ith WT, respectively; ρ is the air density; R denotes the rotor radius; λ i and β i are the tip speed ratio and the pitch angle of the ith WT, respectively; V i is the inlet wind speed of the ith WT; ω r,i is the rotor speed of the ith WT; C p,i and C T,i are the power coefficient and the thrust coefficient of the ith WT, respectively, which can be acquired via the look-up tables or the curves fitting [37].
Besides, the two-degree-of-freedom control framework [29] for each WT is given in  n g is the gear box ratio; τ e and τ p are the time constants for the generator and the pitch angle servo system, respectively.

Wake effect
When the wind inflow goes through the rotated blades, the wind speed behind the rotor will be reduced.Because of this wake effect, the wind speed input of the downstream WTs will be influenced by the operating states of the upstream WTs, as shown in Fig 3 .To describe this influence, the basic Jensen model is adopted to evaluate the downstream wind speed, as follows [38]: where V in is the inlet wind speed of the upstream WT; V x is the downstream wind speed with the distance x directly behind the upstream WT; C T is the thrust coefficient of the upstream WT; D is the rotor diameter; and k is the wake decay constant.
In general, a downstream WT is easily influenced by multiple upstream WTs.By considering the shadowed area by the upstream WTs, the inlet wind speed of the ith WT can be calculated based on the Jensen's multiple wake model, as [39] V i ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where V 0 is the natural wind speed; V ij denotes the wake wind speed of the jth upstream WT at the ith WT location; α ij denotes the intersection weight of the shadowed area at the ith WT location by the jth upstream WT; D ij denotes the distance along the wind direction between the ith and jth WTs; and S ij is the the shadowed area at the ith WT location by the jth upstream WT.Note that D ij and S ij are associated with the dynamic wind direction and the fixed locations of WTs, which can be updated through the coordinate rotation [40] and the overlap evaluation [41] under different wind directions.

Group division of WTs
In this work, the WTs with wake interaction are clustered into a group, while different groups can be decoupled without wake interactions.Based on the wake effect described in Eqs ( 5)-( 7), a large intersection weight α ij indicates a strong wake effect from the jth upstream WT to the ith WT.Hence, two WTs can be put in the same group when their intersection weight is large than a certain value.To simplify this process, the WTs connected in the same transformer are divided into the same group [42] due to their short geographical distances and large shadowed areas, as illustrated in Fig 3.

Optimization model 2.5.1 Objective function.
In response to the AGC signal, the primary task of HPC is to trace the AGC signal as accurately and quickly as possible.Firstly, the power output of all the WTs should be equal to the AGC signal.Secondly, the variations of pitch angles should be as small as possible due to the slow response speed of RSC.Therefore, the objective function of HPC consists of two components, including the power deviation and pitch angle variation, as follows: where f is the overall objective function; μ 1 and μ 2 are weight coefficients for the power deviation and the pitch angle variation, respectively; DP D WF denotes the power variation of the whole WF; Db D WF denotes the pitch angle variation of the whole WF; P ref WF is the AGC signal of the whole WF; k denotes the kth control interval of AGC; and n is the number of WTs in a WF.
2.5.2Constraints.Except the constraints of wake effect in Eqs ( 4)-( 7), the constraints of HPC mainly include the power capacity constraint, operating constraints of pitch angle and rotor speed, as follows [43,44]: where P rate is the rated power of each WT; β min and β max are the minimum and maximum pitch angles, respectively; db dt À � min and db dt À � max are the minimum and maximum ramp rates of the pitch angle, respectively; ω min and ω max are the minimum and maximum rotor speeds, respectively.

Optimization principle of DDSO
In this work, the presented DDSO consists of two main operations, including the surrogate construction and minimum search [45], as shown in Fig 4 .Firstly, the surrogate can be constructed by a radial basis function (RBF) interpolator based on the sample points.Secondly, multiple new solutions can be directly generated and evaluated through the local search based on current surrogate, in which the objective value of the best sample point will be evaluated and used for updated the surrogate.Therefore, a high-quality optimum can be acquired after multiple iterations with these two operations, while a high computation cost can be effectively avoided since the surrogate can accurately replace the original objective function.
3.1.1Surrogate construction.At the initial stage, the sample points are generated by a quasirandom sequence because it can quickly cover the optimization space more evenly than the pure random points.To satisfy the operating constraints in Eqs (10)-( 13), the random sample points can be scaled and shifted to limited within the lower and upper bounds.Here, a frequently-used quasirandom sequence called Halton sequence is adopted to generate the initial sample points, as follows [46]: where x p is the pth sample point; F b d ðpÞ represents a van der Corput sequence with a prime power base b d for the dth dimension; d is the number of dimensions; and n 0 is the number of initial random sample points.
Based on the generation solutions in Eq ( 14), the objective function of each solution can be evaluated according to Eqs ( 8) and ( 9) with the wake effect.Given these data points {(x p , f p )| p = 1,2,� � �,n 0 }, a RBF interpolator is employed to construct the surrogate due to its high flexibility and computation efficiency.Hence, the surrogate can be described as follows [47]: where λ p denote the RBF coefficient for the pth point; k�k 2 is the Euclidean norm; ϕ represents the radial function; and ρ(x) denotes the polynomial function for x.Here, a cubic radial function and a linear polynomial function are selected to construct the surrogate in Eq (15), as follows [48]: where c2R d denotes the coefficient vector of the linear term and c 0 is the constant term.
With the sample data points, the unknow coefficients in Eq (15) can be determined through solving the linear equations, as follows [47]: where λ ¼ ½l 1 ; l 2 ; � � � ; l n 0 � T is the RBF coefficient vector; c p = [c T c 0 ] T is the coefficient vector of the linear polynomial function; ψ is a n×n matrix with T is the objective function vector.
Note that the surrogate can be updated with the added sample points, while the unknow coefficients in Eq (18) will be unique if the inverse matrix is nonsingular [47].

Minimum search.
To achieve an efficient optimization, a local search based on the best point (i.e., incumbent point) found so far since the last reset of surrogate is implemented, as follows [48]: where x q new denotes the value of the qth dimension for the new candidate solution, q ¼ 1; 2; � � � ; d; x q best denotes the value of the qth dimension for the incumbent point; r nd is a random value following the normal distribution N ð0; 1Þ; ε is a random value uniformly distributed from 0 to 1; ε 0 denotes the perturbation probability; x q up and x q low are the upper and the lower bounds of the qth dimension, respectively; and Δδ is the search scale, which will be doubled after three successful searches and be halved after d unsuccessful searches since the last update.
Note that each candidate solution in Eq (22) should be limited within its lower and upper bounds after a perturbation based on the incumbent point.In this work, the number of candidate solutions is set to be 500d at each iteration.To avoid trapping into a low-quality optimum, a merit function is used to evaluate the candidate solutions instead of only using the surrogate.Particularly, it is a linear weighted sum of the scaled surrogate and the scaled distance, as follows [49]: where μ is the weight coefficients of scaled surrogate, 0�μ�1; S(x) denotes the scaled surrogate; L(x) denotes the scaled distance; s max and s min are the maximum and minimum of surrogate among all the candidate solutions, respectively; l max and l min are the maximum and minimum distances among all the candidate solutions, respectively; and l(x) denotes the distance function between different candidate solutions.As a result, a candidate solution with the smallest merit function will be chosen as the sample point at the next iteration, i.e., x best can ¼ argmin p¼1;2;���;500d f merit ðx p can Þ.It will be evaluated by the original objective function f ðx best can Þ, and then used to update the surrogate.As a whole, the execution procedure of DDSO is illustrated in Algorithm 1, where k max is the maximum iteration number; l in max is the maximum distance between all the current candidate solutions and the incumbent point; and l reset is the tolerance distance to stop the minimum search and reset the surrogate.

Power control among different groups.
Firstly, the AGC signal of the whole WF will be distributed to different groups in proportion of their regulation capacities.The regulation capacity of each group can be updated as the wind speed or direction varies, which can be determined via a maximum power optimization, as follows: s:t: Eqs ð4Þ À ð7Þ and Eqs ð10Þ À ð13Þ ð27Þ where n m is the number of WTs in the mth group; P m is the total mechanical power of the mth group; and P i m is the mechanical power of the ith WT in the mth group.The maximum power optimization shown in Eqs ( 26) and ( 27) for each group is solved by DDSO, which can be executed offline with the updated wind speed and direction.Consequently, the AGC signal of each group can be written as > > > : where P max m is the current maximum power of the mth group; ΔP(k) denotes the regulation power command of the whole WF at the kth control interval; P up m and P down m denote the current regulation-up and regulation-down capacities of the mth group, respectively; and M is the number of WT groups.

Power control among different WTs in each group.
To match the objective function of the whole WF, the objective function of each group can be described as follows: ( where f m is the objective function of the mth group; DP D m denotes the power variation of the mth group; Db D m denotes the pitch angle variation of the mth group; and b i m is the pitch angle of the ith WT in the mth group.
The optimization in each group is also solved by DDSO, but it is executed in real-time under the dynamic AGC signals.

Execution procedure.
Taken together, the execution procedure of DDSO based HPC for a large-scale wind farm is provided in Fig 5, where the detailed execution procedure of DDSO can be found in Algorithm 1.The main novelty of this work is included in the offline optimization for regulation capacity, the power control among different groups, and the power control among different WTs, as highlighted in Fig 5.

Case studies
In the case studies, a large-scale WF with 100 WTs is designed to evaluate the performance of the proposed technique, where the geographical distribution of all the WTs is given in Fig 6 .Besides, the parameters of each WT are provided in Table 1.The initial pitch angle and rotor speed for each WT are set to be 0˚and 3.2 rad/s, respectively.The weight coefficients μ 1 and μ 2 in Eq (8) are set to be 1000 and 0.001 as the power deviation is more important than pitch angle variation for AGC signal response.In order to further verify the proposed DDSO, two meta-heuristic algorithms including GA [32] and PSO [33] are adopted for performance comparison, where the main parameters of these algorithms are given in Table 2.

Study of constant wind speed
In this case study, the input wind speed is set to be 12 m/s with the direction θ = 0˚.Moreover, AGC signal for the whole WF is divided into a step signal and multiple continuous step signals.The optimization process of searching the maximum power of each WT group is given in Here, a step AGC signal for the whole WF is set to be 50 MW at the 10 th second.Similarly, DDSO can rapidly find a high-quality solution for each group after the random sample generation, as illustrated in Figs 8 and 9 provides the response process of wind farm for a step AGC signal with a constant wind speed, where DDSO-I denotes the DDSO without considering the pitch angle variation in Eq (8); and

Study of dynamic wind speed
In this case study, the AGC signals are set to be the same as that in the case study of constant wind speed.In the meantime, the wind speed and the wind direction will vary during the power regulation for AGC.

Performance under a step AGC signal.
Here, the wind speed is set to be 12 m/s, while the wind direction changes from 0˚to 10˚at the 25 th second.Although the AGC signal of the whole WF is constant after the 10 th second, the AGC signals to all the groups can adapt

Study of actual power deviation
Table 3 provides the actual power deviation obtained under different wind speeds with different AGC signals, where the actual power deviation can be calculated as P N k¼1 jP ref WF ðkÞ À P act WF ðkÞj; P act WF ðkÞ is the actual power regulation output of the whole WF at the kth control interval; and N is the number of the control intervals for each case studies.Firstly, it can be found that the actual power deviation obtained by DDSO-I is the largest at each case study, which verifies that the lack of pitch angle variation will result in a large actual power deviation.Secondly, DDSO-II can acquire the smallest actual power deviation, which is only 70% of that obtained by DDSO-II under the constant wind speed with a step AGC signal.It also demonstrates that DDSO-II can acquire the higher quality dispatch schemes than that of GA and PSO within a limited number of objective function evaluations.

Conclusions
In summary, the main contributions of this work can be written as follows: 1.The proposed HPC can effectively reduce the optimization complexity and difficulty of power dispatch for a large-scale WF.The consideration of pitch angle variation in the objective function of HPC can make the whole WF rapidly trace the assigned AGC signal.
2. Compared with the meta-heuristic algorithms (e.g., GA), the presented DDSO can acquire a higher quality dispatch scheme for HPC within a limited number of objective function evaluations.Besides, the minimum search can guarantee an efficient local search for DDSO.
3. Four case studies under different wind speeds with different AGC signals are undertaken to verify the proposed technique.The simulation results show that the proposed technique not only can achieve an accurate and rapid trace for different AGC signals of a large-scale WF, but also can adapt to the wind speed variation by fleetly regulating the dispatch schemes.For all the case studies, the proposed DDSO-II by taking the pitch angle variation into account can acquire the smaller actual power deviations than that of GA, PSO, and DDSO-I for the whole WF.Under the constant wind speed with a step AGC signal, the actual power deviation decrement is up to 30% by DDSO-II against to DDSO-I.

Fig 2 .
The control inputs are assigned from the control center of WF, which includes reference active power P ref , reference rotor speed o ref r , and reference pitch angle β ref .The control outputs contain the electromagnetic torque of generator T e and the pitch angle β.The control process is described by a first-order inertia transfer function with a bound limiter and a ramp constraint, where T ref g and T ref e are the reference generator torque and electromagnetic torque, respectively;

Fig 7 . 4 . 1 . 1
It is obvious that DDSO can implement an efficient local search after the first surrogate construction based on the random samples.The obtained maximum power capacities of Groups #1 to #5 are 16.61, 13.01, 11.09, 12.05, and 13.94 MW, respectively.Performance under a step AGC signal.

Fig 9 .
Fig 9. Response process of wind farm for a step AGC signal with a constant wind speed.(a) AGC signal of each group, (b) Power regulation of each group obtained by DDSO, (c) Power regulation of each WT in Group #3 obtained by DDSO, and (d) Power regulation of the whole WF.https://doi.org/10.1371/journal.pone.0291383.g009

Fig 10 .
Fig 10.Response process of wind farm for multiple step AGC signals with a constant wind speed.(a) AGC signal of each group, (b) Power regulation of each group obtained by DDSO, (c) Power regulation of each WT in Group #3 obtained by DDSO, and (d) Power regulaiton of the whole WF.https://doi.org/10.1371/journal.pone.0291383.g010

Fig 11 .
Fig 11.Response process of wind farm for a step AGC signal with dynamic wind speeds.(a) AGC signal of each group, (b) Power regulation of each group obtained by DDSO, (c) Power regulation of each WT in Group #3 obtained by DDSO, and (d) Power regulaiton of the whole WF.https://doi.org/10.1371/journal.pone.0291383.g011